Single neuron computation: from dynamical system to feature detector 
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Abstract 

White noise methods are a powerful tool for characterizing the computation performed by neural systems. 
These methods allow one to identify the feature or features that a neural system extracts from a complex 
input, and to determine how these features are combined to drive the system's spiking response. These meth- 
ods have also been applied to characterize the input/output relations of single neurons driven by synaptic 
inputs, simulated by direct current injection. To interpret the results of white noise analysis of single neurons, 
we would like to understand how the obtained feature space of a single neuron maps onto the biophysical 
properties of the membrane, in particular the dynamics of ion channels. Here, through analysis of a simple 
dynamical model neuron, we draw explicit connections between the output of a white noise analysis and the 
underlying dynamical system. We find that under certain assumptions, the form of the relevant features is 
well defined by the parameters of the dynamical system. Further, we show that under some conditions, the 
feature space is spanned by the spike-triggered average and its successive order time derivatives. 
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1 Introduction 

A primary goal of sensory neurophysiology is to understand how stimuli from the natural environment are 
encoded in the spiking output of neurons. A useful tool for performing this characterization is white noise 
analysis Marmarelis and Marmarelis, 1974 |Wiener, 1958[ Rieke et al., 1997] , whereby a system is stimu- 
lated with a randomly varying broadband signal and the relevant features driving the system are determined 
by correlating output spikes with the signal. Extensions of this analysis to second order have been used to 
find multiple stimulus features to which the system is sensitive |de Ruyter van Steveninck and Bialek, 1988[ 



Bren ner et al., 2000} |Bialek and de Ruyter van Steveninck, 2005 . Recently, these methods have been ap- 



plied to characterize the computation performed by a single neuron on its current inputs | Agiiera y Areas et al., 2 003 , 
Agiiera y Areas and Fairhall, 20"03]|Slee et al., 2005] . 



Due to detailed knowledge of the dynamics of ion channels, it is possible to build dynamical models 
of current flow in neurons that accurately reproduce the voltage response of single neurons to current (or 
conductance) inputs. White noise analysis provides us with the capacity to reduce this complex set of 
dynamical equations to a simple functional model with concrete components: linear feature selection followed 
by a nonlinear decision function generating spikes. This procedure is compelling as it provides a clear intuition 
about the form of the changes in the stimulus to which the system is sensitive, and how the system combines 
these features in the decision to fire. White noise methods have given insight into coding properties at the 
systems level, in visual |Brenner et al., 2000 Fairh alTet al., 2006] |Touryan et al., 2 002, R ust et al., 2005] 



|Horwitz et al., 2005] and somatosensory Petersen, 2003 Maravall et al., 2006 cortex. More recently, these 
methods have also been applied to characterize neural coding in single auditory Slee et al., 2005] , central 
|Powers et al., 2005 and model |Agiiera y Areas et al., 20"03 Agiiera y Areas and Fairhall, 2003] neurons. 



However, mostly missing from these analyses is a clear link between the obtained stimulus features and 
the biophysical properties of the circuit or neuron. This issue is particularly well-posed for single neurons 
under stimulation by direct somatic current injection, where the biophysical parameters of the soma must 
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dominate the neuron's computation. Thus, given the power of white noise methods, the questions that we 
wish to address here are: how does the dynamical system governing neural behavior map to the neuron's 
functional characterization derived using white noise analysis? How are the features determined by the 
neuron's biophysical properties? 

1.1 Minimal spiking models 

Neural dynamics are typically described by conductance-based models which describe the temporal evolution 
of the voltage V due to an input current I and the variable ionic conductances of ion channels: 

CdV/dt = -Y J 9i[V){V -Et)+I, (1) 

i 

where C is the membrane capacitance. Each conductance type i has a reversal potential Ei, where the 
conductances gt corresponding to different ion channels may be voltage dependent through the dynamics of 
the corresponding activation and inactivation gating variables & and fa : 

<?» = Ul o; • (2) 
dti/dt = fafoV) (3) 
dfa/dt = f^ifaV), (4) 

where gi are constants, and the functions and are affine in £ and <p respectively, but not in V . 
The exponents pi and qi are usually taken to be integers, and describe the cooperativity of molecules 
required to gate the ion channel. This set of equations is highly nonlinear due to this cooperativity and 
the voltage dependence of the conductances. While ultimately we would like to understand systems of 
many dynamical variables, we will begin by analyzing a simplified spiking system, the FitzHugh-Nagumo 
(FN) neural model |Fitzhugh, 1 961, Nag umo et al., 1962| . The FitzHugh-Nagumo system is cubic in V, the 
lowest order nonlinearity supporting spike-like behavior, and is defined by: 

ipdV/dt = V(l — V)(a + V) —W + c + 1, (5) 
dW/dt = V-bW, (6) 

where V is voltage, W is a generic inactivation variable with linear recovery dynamics, I is the input 
current and a, b, c and ip are parameters, tp <C 1. The nullclincs of the system are the curves along which 
dV/dt = (the V nullcline, a cubic) and dW/dt — (the W nullclinc, a straight line). Fig[T^ shows the 
nullclines for zero input / on the phase plane. Trajectories have been traced forward and backward in time 
through a randomly chosen set of initial conditions. The unique intersection of the nullclines is the stable 
fixed point. Elsewhere in the plane, trajectories always spiral counterclockwise, but qualitatively there are 
two kinds of orbits, ultimately both converging to the stable fixed point. Subthreshold orbits make only 
small excursions; however, for trajectories starting either with low inactivation (roughly W < —.45) or high 
voltage (for example, V > 0), the system makes a large excursion around the right branch of the V nullcline, 
corresponding to a spike, before settling back to the fixed point. Although we will generally be considering a 
time- varying input I(t), we will make heavy use of the 1 = phase plane: one can regard this as representing 
the instantaneous flow if the input current is switched off. For 1 = 0, one can define a separation of spiking 
initial conditions from subthreshold conditions, marked in light vs. dark gray in Fig[l}i. This separation 
of behavior can be thought of as a threshold. Typically the threshold for spiking is taken to be a threshold 
in voltage only, but this phase plane picture demonstrates that one should really consider the threshold to 
depend on both V and the hidden inactivation variable(s). 

One can approximate this "dynamical" threshold by considering a minimal spiking trajectory. There is no 
absolute distinction between spiking and subspiking trajectories in type II neurons |Gerstner and Kistler, 2 002] 
such as Hodgkin-Huxley and FitzHugh-Nagumo, so this choice is necessarily somewhat arbitrary. Starting 
with a point on the phase plane at the end of the spike plateau, one can eliminate time from the dy- 
namical equations, solving parametrically for the curve passing through this point. Here, we obtained one 
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such curve by numerically solving the equations from the local maximum of the cubic nullcline, given by 
^V(l-V)(a + V) = 0. _ 

One of the issues we will explore in this paper is the relationship between the geometry of the threshold 
in the multidimensional space of the dynamical variables and the measurements one makes in white noise 
analysis. 



1.2 White noise analysis 

Our approach will be based on the functional perturbation expansion method, which we will briefly review 
here. The simplest such expansion is known as the Volterra series |Volterra, 1930 Marmarelis and Marmarclis, 1974 



Marmarelis, 2004] , where the functional y(t) = y[I{t)] is expanded as 



y(t) = ho + Jdtih 1 (t 1 )I(t-t 1 ) + j l J dt 1 dt 2 h 2 (t u t 2 )I(t-t 1 )I(t-t 2 ) + --- , (7) 

where the kernels hi are called Volterra kernels. One problem with Volterra analysis is that the successive 
order kernels are not independent. For example, h 2 contains a constant component ~ J dt\ h 2 (ti,ti). 

Wiener analysis Wiener, 1958, Rie ke et al., 1997] avoids this issue. In this framework, I(t) is assumed 
to be zero mean Gaussian white noise satisfying 

(I(t)I(t'))=a 2 5(t-t'). 

Under Wiener analysis, y[I(t)] is expanded as 



y(t) = g + J dhg^Iit-h) 

+ h{J dt ^92( t ut 2 )I(t~t 1 )I(t-t 2 )-<7 2 J dt' g 2 (t',t')^ + ••• . (8) 

In the Wiener expansion, each kernel is independent from the others, since at each order, one subtracts out 
the lower order components. Therefore, the kernels g% can be easily obtained by measuring the correlation 
functions between y(t) and I(t). For example, 



+^-\ I dt 1 dt 2 g2(ti,t2){I(t-t 1 )I(t-h)}-a 2 I dt' g 2 {t',t')\ + 



(»(*)) = 9o + J *i5i(ti)(J(*-ti)) 
1 

"2! 

= 3o, 

(y(t)I(t-ti)) = gi(h)/a 2 , 
(y(t)I(t-h)I{t-t 2 )) = g 2 (t 1: t 2 )/a\ 

and so on. The correlation method can be used to determine the relationship between Wiener and Volterra 
kernels. From Eq ([8]), 



5o = h+ / dtifci(t)(I(i-ti)) 



+ / </M/2/.'2(/|../2.) ;/-/ - /.) )/(/ - 1-2)) - 

gi(t) = /',•:/; • y I dt'h^l'.l'.i) + 
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and so on. 

Wiener analysis has been successfully used to determine the relevant feature space for neural systems by 
reverse correlation of the spiking output with white noise input. For this application, the output is taken to 
be the sequence of spike times, p(t) = ^\ S(t — ti). Let us define the stimulus s as a linear transformation 
of the current input I(t): 

8i= [ dt' fi(t')I(t-t') (9) 



where fi are some set of independent linear filters. The goal is to determine how a single spike is generated 
by /, reduced to a vector s comprising projections along the relevant dimensions 

We need to compute the probability of spiking as a function of the input, P (spike|s): 

P (spikc|s) = P(spike)G(s), (10) 

where G(s) is the input/output function relating p(t) to the time- varying stimulus s(t), and we have separated 
out the scale factor P(spike) proportional to the mean firing rate. In principle, since the equations of motion 
are deterministic, P (spike|s) is a Boolean function. However, as we will be approximating s by a finite 
dimensional truncation, and will be solving for P itself in successive orders, the result will no longer be 
Boolean. 

We can compute P (spike|s) using Bayes' theorem |Rieke et al., 1997] , 

= P(spike|s) = P(s|spike) 
v ; P(spike) P(s) v ; 

This shows that the Wiener kernels of G(s) are the moments of P(s|spike), which can be easily measured. 
In this paper, we will focus on the first two moment^] of this distribution. The first moment is the spike- 
triggered average (STA), s(t): 

s(t) = (I{U - t)) l = JdssP (s|spike) = / ds sP(s)G(s), (12) 

where the average over the spike occurrence times {ti} is replaced with an average over the spike- 
triggering ensemble. When the neuron has multiple characteristic features, the STA may contain little 
information about them since it only points in a single direction which is a certain linear combination of 
relevant features. The second order kernel, however, can give such information. This kernel is related to the 
spike-triggered covariance (STC) matrix G, which we will define as 

G(i,0 = (I(t t -t)~s(t),I(U-t')-s(t')) l -(I(T-t)I(r-t / )) T 

= a 2 g 6(t,t') + a A g 2 (t,t')-s(t)Kt')-<r 2 t(t,t') (13) 
= a%(t,t') - s(t)Ht')- 

g2(t,t') is the second order Wiener kernel of G(s). When a functional F[s] depends on Gaussian noise s, 
(sjP[s]) = J2j ( s i s j) (djF[ s ])- Thus, we can expand g 2 |Bialek and de Ruyter van Steveninck, 2005| : 

where Eq. ([5]) introduces the linear filters which define the components of s. 

Therefore, unless the modified Hessian — / g s d 9s . G^ — {^^-G^ {^^-G^) has null direction^], the 

STC is spanned by the relevant feature space {/j}. Diagonalizing the STC and extracting the eigenvec- 
tors with largest absolute eigenvalue will determine the leading dimensions spanning the feature space 



1 Note that the zeroth order kernel is simply go = 1- 

2 Since s(t)s(t') is only rank one, there are only two possibilities leading to null directions. The first is when (didjG) has null 
directions, in which case it is excluded (for our purposes) by definition. The second case is when (didjG) has an eigenvector 
s(t) with an eigenvalue ||5(t)|| 2 , which does not happen generically. 
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Brenner et al., 2000]|Aguera y Areas et al., 2003 Bialek and de Ruyter van Steveninck, 2005] . Eigenmodes 



appearing with negative eigenvalue correspond to stimulus directions in which the variance of the spike- 
conditional distribution is reduced relative to the stimulus prior distribution; eigenmodes with positive 
eigenvalues are directions with increased variance. Such cases may arise if the spike-conditional distribution 
is bimodal in some direction, or if the spike-conditional distribution forms a ring around the origin; phase- 
invariant (complex) cells are an example Simoncclli et al., 2004 Rust et al., 2005, Fairhall et al., 2006, Peter sen, 2003| 
Maravall et al., 2006] . 

In this paper, we work with a multidimensional space of dynamical variables, as in Eq. (1-6), for which 
we obtain the Volterra series. When we denote them with j/j, the Wiener kernels for G(y;) can be expanded 
as 



9i[G] 



>^ [Vi. 



QQ d 2 G 

92 [G] = E [Vi] + 22 a„,. a „.. ffi \Vi\9i [Vj] 



dyidyj ■ 



where g n [yi] is a nth Wiener kernel of and dots represent higher order terms than quadratic. Here we 
note that g n [yi] is generated by the infinite series of Volterra kernels. Moreover, we will see later that the 
Volterra kernels at each order can be constructed from the set of first order kernels. 



2 Linear systems 

In the previous section, we discussed the application of white noise methods to experimental data where 
the output of the system is measured as the firing times of spikes. With access to the entire dynamical 
system, one could attempt to model the relationship between the input I(t) and the output voltage V(t), 
including both subthreshold and spiking behavior. As we have discussed, neurons are highly nonlinear. 
Perturbative expansions such as Wiener/ Volterra series generally are unable to capture the nonlinearities 
of spike generation with few terms. However, a very successful and highly influential approach to the 
modeling of neural systems has been to separate a linear filtering stage from an explicit nonlinearity cap- 
turing the very sharp voltage increase of the spike [Victor and Shapley, 1980] |Gerstner and Kistler, 2002"] 
|Kistler et al., 1997] |Berry II and Mciste r7 1999] . This is known as a cascade model, and in various forms is 
the theoretical basis for many of the simple neuron models in use, including the integrate-and-fire neuron 
Gerstner and Kistler, 2002, Agiiera y Areas and Fairhall, 2003 , the Spike Response Model Gerstner and Kistler, 200*2] , 



generalized integrate-and-fire models |Paninski et al., 2004| and cascade models of retinal ganglion cells 



|Victor and Shapley, 1980] |Victor and Shapl ey, 1979] |Keat et al.,"200l| |Pillow et al., 2005| . The success of 



these models suggests that a breakdown of the dynamical system into its linear component and a threshold 
is a useful approximation to examine. In our analysis of the dynamical system, we will take the approach of 
experimental work, where the output is reduced to a sequence of spike times. Our goal here is to derive the 
form of the multidimensional linear stage that arises directly from the equations of motion, and to examine 
the consequences of thresholds of various forms. 

We will begin by walking through some simple linearized systems analysis and make clear the connection 
with the cascade model picture. Rewriting Eqs. (1-4) in the following simple form, 



dtVk = fk(yi,V2, ■ ■ ■ ,Vn) + I(t)S k0 , k = 0, l,--- ,n- 1, (15) 
of the system around the fixed point is given bjO 



we assume that they have a fixed-point solution?^ (t) = j/i when I(t) = 0. Now the linear approximation 



V km y^ = I(t)S k0 , (16) 



5 From here on, we use the Einstein summation convention, e. g. Cfc m <i m = J^m Cfcm^n 
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where 

g f 

T^km — &km &t J km-, J km 



(17) 

(0) 



dy 

The linear filters are the kernels of Eq dTTJ) , e (1) = 2? _1 . They satisfy 

V kr J^(t) = 5 kl S(t), (18) 
giving the linear filter form for the system's evolution: 

y£\t) = f dt'e^it-t'W). (19) 



where e^ 1 ' = e^ 3 . Eq (JIBJ can be solved by diagonalizing the Jacobian matrix J km . Let's consider the 
following diagonalization with eigenvalues X a , 

A Q <5 Q /3 = U a kJkmU ni i. (20) 



Here we employ the notation that indices in the diagonalized basis are Greek letters. Therefore, Eq (|T8|) can 
be written in the new basis as 

Vj£ = U a0 S(t), V a = d t - X a , (21) 
where the first order kernels e^' in this basis are 

eU(t) = e x * t H(t), 

where H(t) is the Heaviside step function and arises from causality. The eigenvalues A a can be real or 
complex. The sign of the real part of the eigenvalues is determined by the stability properties of the fixed 
point at j/ ' . For a neural system this fixed point should be attracting, so all eigenvalues have negative real 
part and thus eigenvectors decay as t — > oo. This underscores that there are only a limited number of forms 
that the linear filters can take: decaying exponentials or damped oscillations. 

We note the following facts. First of all, in general, an n-dimensional system will have n independent 
first order kernels, except when the Jacobian has a degenerate spectrum. Further, the linear kernels can be 
generated from a "master kernel" 

a 

and its time derivatives up to nth order. To show this, from the expression of the derivative of a master 
kernel 

dt k tp(t) = x t eXat = T ka e x "\ T ka = A*, (k, a = 0, . . . n - 1), 

a 

up to additive singular terms coming from the derivatives of H(t). T ka is a Wronskian matrix of exponentials 
at t = 0. The determinant of T ka is the Vandermonde determinant det(Tfco,) = ri Q >/3(Act — A/3), which cannot 
vanish by definition Gradshtcy n and Ryzhik, 2000] . This means that the master kernel and its derivatives 
can be transformed to e Aat 's in a non-singular way. 

Therefore, we conclude that all first order kernels can be expressed in terms of a master kernel and its 
time derivatives. Frequently the first order kernel is non-zero at t = 0, so that the time derivative of the 
filter also includes a delta function at t = 0. We will see the appearance of this singular component in the 
white noise analysis. 

To recover the filters corresponding to the dynamical variables, one must invert the diagonalizing opera- 
tion, Eq. P0|) to obtain linear combinations of the eigenmodes: 

4 1] (t) = U- a ^\t) =J2c a e^H(t), (22) 



G 



where the coefficients c a derive from the components of the Jacobian. Thus the filters of the linear system are 
sums of exponentials, either purely real or with an imaginary component. To the extent that subthreshold 
dynamics are well approximated by the linearized system, this shows clearly why one might expect to find 
both integrate-and-fire-like neurons, with filters having purely real associated eigenvalues, and oscillate-and- 
fire neurons [izhekevic h, 2001] |Hutcheon and Yarom, 2000] with a corresponding eigenvalue having nonzero 
imaginary part. In terms of the linear system, the possible set of filters or feature space of the system is dual 
with the dynamical variables which define the system. The subset of features that are relevant to a spike 
occurring are those that contribute to crossing the internal state over "threshold" . 

The properties outlined above have an interesting implication for white noise analysis in the case that the 
system is well described by subthreshold linearity. As we will see in the next section, the STA is typically 
a linear combination over the relevant features and their derivatives. If the STA has components onto all 
the true features, as it typically will, one can derive the basis for the feature space from the STA and its 



successive order time derivatives, as has been empirically observed Rieke, 1991 



2.1 Higher order series terms 

The higher order approximations to the system can also be expressed in terms of first order kernels. For 
example, the second order approximation is given by the following equation: 



T> ( 2 ) — — H W f 1 ) n — d 2 /fc 



(23) 



whose solution is 



(2) 

Vk 



^ J ds 1 ds 2 £ { k\t- si,t- s 2 )I(si)I(s 2 ), 

e k 2) (t x ,t 2 ) = J dt' H lmn e^(t')e^(t x -t')e^(t 2 -t'). (24) 

In this way, all of the higher order kernels can be obtained by outer products of first order kernels. The 
detailed computation is summarized in appendix [5] 

2.2 White noise analysis of threshold-crossing linear neurons 

In this section, we will assume that the higher order nonlinearity of the system is well captured by a threshold. 
Since voltage is generally the only observed variable, "the threshold" is often taken to be a threshold on the 
voltage. The imposition of a threshold on the voltage alone has been applied even to neural models such 
as FN and the Hodgkin-Huxley model system, where there is no difficulty in defining a threshold in V and 
the "hidden" dynamical variables of W (FN) or to, n, and h (HH). Here we will consider the implications 
of considering the threshold to apply in multiple dimensions. Before we discuss the validity of the threshold 
approximation for a FN neuron, we will illustrate the results from applying white noise analysis to models 
with a variety of threshold structures acting on subthreshold linear dynamics given by the linearized FN 
system. 

Previous work has treated simple one-dimensional linear/nonlinear models, in which a spike is defined 
when the output of a single linear filter / on a random current input crosses a threshold. One can show 
DcWccse , 1995| [Meister, M., private communication] that covariance analysis on this model finds two modes, 
the filter / and its time derivative, /'. The derivative mode is a result of the criterion defining the spike only 
on the upward crossing of the threshold: 

d ft r* 

drf(t - t)I{t) > -» / dTf'(t-T)I(r)>0. (25) 



dt 

Thus the variance of projections onto /' is also reduced with respect to the prior. Covariance analysis has 
also been performed |Fairhall et al., 2006] on more realistic models in which an exponentially decaying af- 
terhyperpolarization is added to the voltage following a spike Gerstner and Kistler, 2002, Kc at et al., 2001] . 
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This makes the threshold-crossing model more realistic, but the manipulation does not alter the (/, /') eigen- 
mode structure unless the timescale of the AHP is short enough that spiking occurs many times during the 
superthreshold fluctuation, which destroys the correlation with positive /'. 

Here we will assume that the system evolves according to linear subthreshold dynamics. We treat five 
cases. In the first three, the threshold is linear, but in different components of the 2D space; in the last two, 
the threshold is a nontrivial 2-dimensional structure. The thresholds used are: 

1. a traditional threshold in V only 

2. a threshold in W only 

3. a threshold on a linear combination of V and W 

4. the union of two piecewise linear segments in V and W , i.e. a logical 'or' on I and 2 

5. a curved threshold on a smooth, nonlinear function of V and W . 

From the discussion in <J21 the linear space in which the n-dimensional system operates is defined by the n 
filter inversions of the dynamical system, Eq. (fH))) . The effect of the threshold is to select as relevant from 
those n dimensions the subspace defined by the filter directions that have a nonzero projection onto the 
threshold. Thus we would predict that for the linear threshold, we will find two relevant modes, the primary 
direction / which is normal to the direction of the threshold, and /', its time derivative. For a threshold 
in a single dimension, V or W only, / should recover the V or W filter respectively, and the corresponding 
time derivative. For a threshold that is a linear combination of V and W, f should be a linear combination 
of the V and W filters. The distribution of trajectories is simply a linear transformation of the original 
Gaussian stimulus distribution. Any linear threshold in this plane will produce negative eigenmodes, as the 
set of points selected by the threshold must have decreased variance with respect to this prior. 

When the threshold is not linear in V and W, we expect to find two primary filters /i and fa which 
will be some linear combination of the V and W filters, and in principle, the time derivatives. However, the 
number of modes will be less than or equal to n + 1 because, as we have shown, the time derivatives span 
the same space as the linear variables themselves. (Note that the +1 accounts for the singular component of 
the time derivative, 8(t), arising from discontinuity of the filters at t = 0.) Furthermore, a threshold which 
is not linear may produce a spike-conditional distribution with a direction in which the variance is actually 
increased with respect to the prior; such a direction will appear with a positive eigenvalue. 

Fig shows the thresholds that we applied to the FN first order linear dynamics. The results of the 
covariance analysis are shown in Figs[2j; and d. Fig[2j: shows the covariance eigenvalue spectra for all five 
cases. The significant eigenvalues are empty circles beyond the error level denoted by the shaded box0 In 
the cases with thresholds on V only, on W only, and linear in both V and W, only two modes are obtained. 
With 2D thresholds, cases 4 and 5, three eigenmodes appear. Fig[3H shows the corresponding modes. We 
see that the V threshold picks out 6y, the W threshold selects e[^, and the (V, W / )-linear threshold leads 
to a linear combination of these filters. In the latter two cases, the eigenmodes are a linear combination of 
ty\ eyj) and dt€y \ itself a linear combination of the V and W kernels and 5(t). 

Fig [Hi also shows the projection of spike triggered stimuli onto two of the distinguished eigenvectors. In 
the first three cases, the linearity of the threshold is manifested in the 2D plane as a single elliptic cluster, 
while the non-linear threshold cases show richer structure. 

Recall that the filters are generated by linear combinations of a master filter and its time derivatives. This 
implies that our multidimensional threshold structure has a natural interpretation as a dynamical threshold, 

4 Significance of eigenvalues can be evaluated as follows. Let N be the number of spikes and d the dimension of each sam- 
pled stimulus. To estimate the finite size/dimension error, we choose N d-dimensional stimuli at random |Rust et al., 2005 . 
In this case, the covariance matrix C in Eq. I|14| l is filled with d(d + l)/2 random numbers drawn from a normal distribu- 
tion JV(0, cr 2 v/2/iV) in the large N limit. When d> 1, the eigenvalues of this matrix follow Wigner's semicircle law, and their 
distribution is bounded by ±2v / 2<3' 2 y'd/N, which is the error level |Mehta, 1991] . When the stimulus has small correlation, the 
semicircle is not a good approximation of the eigenvalue distribution, but change in the upper and lower bounds rarely exceeds 
an order of magnitude, which can be checked by numerical simulation. 
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depending not only on the voltage V , but also on dtV, d 2 V ', and so on |Azouz and Gray, 2000 . We will 
show how well this approximation fits for the Hodgkin-Huxley model in the final section. 



2.3 Analysis of threshold-crossing models in multiple dimensions 

An advantage of the simple models introduced in the previous section is that they can be treated analytically. 
From Eq. we have 

_ P(spike|s)P(s) _ P(spike|s)P(s) 
^ ~ P(spike) ~~ / P>sP(spike|s)P(s)' 

where J Ds = J ds\ds2 ■ ■ ■ ■ Instead of directly computing G(s), we characterize it by a moment generating 
function, 

W\j] = log J J DsP(spike|s)P(s)e j s . (26) 

We define an n-dimensional first order system, given by dynamical variables y k and the corresponding 
linear filters e k as in Eq. (|16[) . We assume that the stimulus is a Gaussian white noise current I(t) with zero 
mean and variance a 2 . Thus as before, 

roc 

y k (t) = / dre fe (r)/(t-r). 



We denote a random segment of the current stimulus 7(r),r < t as an infinite-dimensionajf) vector s. 
Any sample of y is therefore a functional of the random variable s, y[s]; for simplicity we will simply write 

y- 

Spiking is determined by crossing a threshold 9(y) = in the phase space from below. Then, 

P(spike|s) = S(0(y))H(y ■ V y 0(y))y • V y 0(y), 

The Heaviside function H(-) ensures that spiking only occurs on a threshold crossing from below, and the 
weight factor y • V y 0(y) is a geometric factor accounting for the flux at the threshold. The filters e k are not 
necessarily normalized or orthogonal to each other, and it is convenient to define the orthonormal basis {/^j 
which spans the same space as {e k }- Then there will be a linear transformation T^ k 

and so we define a new coordinate system 

As s is uniform Gaussian, the orthonormally transformed variable z is also uniform Gaussian with variance 
a 2 . 

Now we can separate the stimulus into two components: its projection into the subspace spanned by the 
{ffj,} an d the orthogonal component. We will denote the corresponding directions of j as j = j|| + where 



J'||^(t) = / drj(t-T)/ p (r). 
Jo 

The moment generating function can then be separated as W[j] = W\\ [jn] + Wj_[jj_] where 

/2 
Ds± e sl/2^+i x .s x = 1-j^t) 2 + const, (27) 

5 In practice, any application to data requires discretization in time. We use the convention that a function / (t) is discrctizcd 
as ft = f(t)V At where At is the time step. This implies J f 2 dt S3 f(iAt) 2 At = f 2 = f • f and thus the vector f obtained 
by discretizing f(t) has a vector norm corresponding to the L 2 norm. 
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where s_l are the components of s orthogonal to the plane spanned by z, and 



w \\lh 



log 



log 



/ 
/ 



d n z det(T- 1 )e- z2 / 2CT2 8{6{z))H{z • V z 6»(z))z • V z 6>(z)e j ir z 



d n z e- z2/2a2 5{9(z))w{z)e^ - z + const, 



(28) 



with 



w{z) = H((Rz) ■ VJ(z))(Rz) • V a 0(z), 



where R = TJT 1 is the Jacobian matrix of the linear system, Eq. (fTB|) . in the {f^} basis, and we use 



As previously discussed, the closure under time differentiation of the space spanned by the first order kernels 
is assured by z, as in Eq. pTj) in the appendix. Note that our separation of W[j] depends on this particular 
property. For a fixed spike time t, restricting ourselves to the region where r < t, W± [j] becomes a static 
integral. Thus Eq. (f2"g)) shows clearly how the properties of this model emerge. P (spike|s) is just given by a 
threshold with a weight function w{z) up to some linear transformations. Therefore, for example, the linear 
threshold case reduces to the one dimensional case. Since V z #(z) is a constant vector, every computation 
reduces to a one dimensional integral in this direction and the corresponding filter is a linear combination 
of f^s, as we have observed in Fig Oi. Furthermore, a variety of structures may be generated depending on 
how trajectories cross the threshold. 

From Eqs (|28p and (f2"T)) , we can derive the analytic forms of the first and second order moments: 



Note that we have assumed a prior based on the distribution of randomly driven trajectories of the linear 
system. This takes no account of perturbations of this distribution due to flux from the system's return from 
a spike. This assumption is equivalent to assuming that the last spike is in the distant past, so that memory 
of that perturbation has vanished. In this paper we will consider only this "isolated spike" case. We will 
return to this point in the discussion. 

2.4 Variance dependence 

Through Eqs. (f29|) and (|30|) . this model captures an explicit dependence of the white noise outcome on 
the stimulus statistics, in particular, the variance a 2 . While the subthreshold dynamics are linear, the 
dependence on a is nonlinear due to the threshold shape and the weight function w(z). We discuss two 
examples below. 

We first consider a linear threshold. Through a suitable linear transformation, this case simply reduces 
to a filter-and-fire model with a single filter, say /o, and a fixed threshold z$ = Of. Now the distribution of 
threshold crossing points is constrained by zq = Of and zq > 0. As we mentioned in the previous section, 
Zq lies in the originally defined feature space, and is given by another single filter, which we denote by f\. 
In other words, when the normalized /o is denoted by /g, we can choose fi — f' by a suitable orthogonal 
linear transformation. Thus, our system depends on two filters, which are depicted schematically in Fig [3^,. 

The STA is the centroid of the distribution of threshold crossing points. Eq. (f29|) reproduces a previously 
known result [Meister, M., private communication] [Fair hall et al., 2006] , 



Zfj, = T^kVk = T^kJkiVi = {TJT 1 ) filJ z„. 




(29) 



and 




(30) 





/i(r). 



(31) 
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With high variance a ^> 9f, the STA is dominated by f\. In the feature space, with increasing variance the 
threshold stays the same, but a larger portion of it is crossed by trajectories driven by the larger variance 
ensemble, as can be seen in Fig [3^. 

When the threshold is curved, the a dependence is considerably more complicated. We will consider an 
extreme but analytically tractable version of this case to illustrate the point: let the curved threshold be 
approximated by two linear ones imposed at Of and 9 g in the fo and go directions respectively, as in FigO). In 
this case, one segment of the threshold imposes a dependence on the filters fo and its (normalized) derivative, 
fx, while the other selects go and g\ = go/ \\ go\\. The space of relevant features is still two dimensional, or 
three including the S- function, since go and g\ are linear combinations of fo, /i and possibly also S(t). 

Now the STA of this system is 

STA(r) = cos 2 ip ■ STA/(r) + sin 2 ip ■ STA 9 (r), tp = tan" 1 e ( e /-^)/ 4CT2 , (32) 
where STA e (r) = 6» £ e (r) + -^= ei (r) as in Eq. ([3l]). Note that as in Fig [3b, the STA does not lie on the 

' 7r/2 



threshold. As for a variety of experimental examples such as compl ex cells [Tburyan et al., 2002 , neurons of 



rat barrel cortex [Petersen, 2003| , and some retinal ganglion cells Fairhall et al., 2006 , the spike-triggered 



stimuli are poorly represented by their first order statistics, the STA. Also, the coefficients of STA/ iS depend 
exponentially on 9f i9 and a. Thus, the system shows a nonlinear dependence on the stimulus variance. 

Fig[4]shows an example. The W threshold model does not show any significant change, except sharpening 
of the STA due to the increased component of ew, Fig UJd, and linear broadening of the spike triggered 
stimulus distribution. The piecewise linear threshold case is more dramatic: while the smallest variance does 
not drive the system hard enough to produce a positive eigenmode (hence the one-dimensional distribution 
of projections seen in Fig. [4Jd and c, red), a new significant mode emerges as the variance increases. The 
STA changes beyond sharpening and almost looks like a different model at high variance compared with low 
variance. Each of the significant modes change as more trajectories cross the threshold from the other side 
and the principal axis of the distribution rotates, as seen in Fig[4j:. 



3 Dynamical threshold 

In the previous section, we considered neuron models composed of linear filters derived from the FN neuron, 
and some choices of imposed thresholds. However, when we consider the full FN neuron, the threshold arises 
from the structure of the FN equations ([5])- ([6]). Here we discuss the importance of the threshold identification 
for reverse correlation analysis. 

The problem of the identification of a threshold arises immediately upon attempting a reverse correlation 
analysis; spike times are often defined by the threshold crossing of the voltage. Here also, we can impose an 
arbitrary threshold in V to identify each spike. The STA obtained using this scheme is displayed in Fig [5^. 
It is clear that it cannot be well described by the first order kernels. However, this does not mean breakdown 
of the analysis; rather, it underscores the point that the STA or any other single spike quantities should be 
computed using the "correct" threshold that we take to be the dynamical threshold discussed in SJTJ 

FiglSJa compares the spike triggered stimuli in the fixed V and dynamical threshold cases. While the peaks 
of the stimuli are spread out in time in the fixed V threshold case, for the dynamical threshold the peaks lie in 
a narrow band around the spike time. As we can see from Fig [I] the spiking trajectories cross the dynamical 
threshold before crossing the V threshold, and this timing difference, At, depends on W. Additionally, since 
the system is driven by Gaussian white noise, the variability increases with the timing difference, inducing 
a point spread function on the STA. Hence, each spike triggered stimulus is contaminated by this temporal 
jitter, and estimated filters are distorted. As mentioned in fll.ll the FN neuron does not have a clear-cut 
dynamical threshold; the threshold was chosen to some degree arbitrarily, and this will induce some error in 
the estimation of filters. However, Fig[6^ shows the improvement in the STA computed using the dynamical 
threshold. Figs \5fc and d show the point spread distribution - the distribution of values of At - and how 
it is correlated with W. This situation is similar to that discussed in |Aldworth et al., 2005| , where it was 
noted that temporal or spatial jitter can blur the estimation of filters and receptor fields. There, a blind 
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deconvolution algorithm was used to "dejitter" and realign the spike-triggered stimuli, which dramatically 
sharpened the estimate of the spike-triggering stimulus. In a real system, jitter may indeed be due at least 
partly to noise, but it is also possible that a component of such jitter is deterministic and due to variability 
in the point of dynamical threshold crossing, as in the FN case. Blind deconvolution may then be viewed as 
an empirical approach to recovering a dynamical threshold based on spike times originally recorded using a 
voltage threshold. It is interesting to note that estimated jitter in such a case is correlated with projection 
onto the STA derivative!! 



4 Fully nonlinear systems 

In this section, we discuss the covariance analysis of the full FN system and compare the result with the 
same analysis of the first and second order approximation. 

We identified ~2x 10 6 spikes first using a voltage threshold, and then backtraced each trajectory to the 
point where it crosses the dynamical threshold shown in Fig [TJ The trajectory might cross the threshold 
multiple times before spiking due to the noisy input. In this case, we used the first crossing point after the 
trajectory diverges from that of the second order approximation. This is based on the assumption that the 
second order is a sufficiently good approximation of the system in the subthreshold regime. For both the 
first and second order system, crossing of the dynamical threshold is used to identify spikes. Right after a 
spike, we imposed a post-spike inhibitory period equal to about a "spike width" , as empirically determined 
from the full system. 

The results are shown in Fig [S] and [71 Fig [5^ shows that the STA of the first order approximation is quite 
similar to that of the full system, and the second order STA is even closer. From the covariance analysis, 
we see that the FN neuron, like the HH model |Agiicra y Areas et al., 2003 , has one dominant negative and 



one dominant positive eigenvalue. However, the results from the covariance analysis of the three systems 
differ considerably. The second order and the full system both show a relatively large number of significant 
eigenvalues. This is due to the contribution of the second gi\V\ and [W] (and higher order) kernels in 
Eq (TT5]) . However, these are relatively suppressed compared to the first order modes. Further, the spectrum 
derived from the second order system, Fig[6j3, has no positive significant eigenvalue. 

Fig Eh provides more detailed information. The first order case is just as we have seen previously, with 
three modes which are well described by linear combinations of the linear kernels. In the second order case, 
Vi and V2 are comparable to those of the first order. V3, which is not approximated by the first order kernels, 
must arise from the second order kernels. In this regard, the full system is well matched with the second 
order approximation. There also, V3 is comparable to that in the second order, and vi y 2 are derived from 
the linear kernels. However, the full system also has a positive mode, v+, from the linear kernels, which 
resembles the first order rather than the second. 

The geometry of the spike-triggering stimulus projections is shown in Fig [7^,. The first order, not 
surprisingly, recapitulates the curved threshold case in Fig However, the second order is more like 
the (V, FT)-linear threshold case, while again the full system resembles the curved in Fig A possible 
explanation for this is that each system probes different parts of the dynamical threshold. Fig[7jD marks the 
density of threshold crossing points for each model. In contrast to the first order and the full system, which 
access a large section of the threshold with non-trivial curvature, the second order only probes a small and 
almost linear section. This is the reason for the lack of a positive eigenvalue. As in the toy models in <|2j 
the contributions of the relevant dimensions are determined not only by local information (filters) but also 
by the global structure of a multidimensional or dynamical threshold. 



6 This correlation arises because the STA derivative is the linear approximation to time translation of the STA. Hence, if 
a current history aligns with the STA better under a small time translation or jitter, it will have a projection onto the STA 
derivative proportional to that jitter. In Sect. 2 of this paper it was demonstrated that the STA and its derivatives are, up 
to a linear transformation, the filters associated with the dynamical variables of the linearized system. Therefore there is a 
relationship between blind deconvolution to optimize fit to an STA and estimation of a dynamical threshold based on the 
"hidden" (non-V) dynamical variables. 
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5 Abbott-Kepler model 



In this section, we apply the same analysis to a two-dimensional model which is more nonlinear and more 
realistic than the FN model. 

Abbott and Kepler |Abbott and Kepler, 1990| developed a two dimensional reduction of the Hodgkin- 
Huxley model, based upon the observation that there is a separation of timescales between the faster m and 
the slower n and h variables, m is then replaced with its asymptotic value at the membrane voltage V, while 
n and h are controlled by another voltage variable U. The equations of the Abbott-Kepler (AK) model arc 
of the form 

= f(V,U)+I(t), (33) 

II = • 9( ^' t/) ' (34) 

where f(V, U) and g{V, U) are nonlinear functions in V and U. Their derivation is briefly sketched in 
appendix [C] and we refer to the original paper |Abbott an d Kepl er, 1990] for further detail. The nullcline 
for U is given by g(V, U) = 0, which is satisfied by V = U . The V nullcline, f(V, U) = 0, is more complicated 
and obtained numerically. The two nullclines intersect at a fixed point V = U = — 65m V. 

Fig [8^, shows the phase plane of this model with zero input current. Like Fig[TJ the threshold structure, 
which can be obtained numerically, is visible. Due to the strong nonlincarity, spiking trajectories arc well- 
defined on the phase plane and the threshold has less ambiguity than for the FN model. Again, we try to 
identify the dynamics of the system in the subthreshold regime with the first order approximation. Unlike a 
FN neuron, the Jacobian of this system has complex eigenvalues A± = — 0.2118± i0.4035ms~ 1 , and therefore 
the first order kernels 6y] \j oscillate, Fig ITDk . This is consistent with the oscillatory linearized behavior 
associated with the full Hodgkin-Huxley model near equilibrium. 

Before we carry out covariance analysis on this model, we examine the effect of a dynamical threshold, 
as in section [3] Fig [9^ shows the spike triggering stimuli aligned both with a fixed threshold in V, chosen as 
V = — 40m V to unambiguously select spiking trajectories, and the dynamical threshold. The typical time 
shift is of order < 1ms; the overall STA only suffers from a slight time/phase shift when the V threshold 
is used. However, in the small time scale of < 1ms, there is discrepancy from the dynamical threshold 
case, which is characterized by a large delta-function component at the spike time. Both STAs are well 
approximated by linear combinations of ejj^, up to a small deformation around -15ms due to the effects of 
multiple spikes. 

Fig Hni shows the results from covariance analysis carried out with ~ 2 x 10 6 spikes for Gaussian white 
noise stimuli with various variances. For comparison with previous results, we selected three eigenmodes 
corresponding to the leading two negative and the largest positive eigenvalues. We find that they are reason- 
ably well approximated by linear combinations of £y\j in most cases, although they are sometimes affected 
by a ^-function at t — and a large multi-spike effect at high variances. The multi-spike effect can be elimi- 



nated by considering only the isolated spikes when the spike rate is low Aguera y Areas and Fairhall, 2003 
At higher variances isolated spikes are rare and there is a stronger influence of oscillating "silence modes" 
|Aguera y Areas et al., 20"03 . As in the FN neuron case, we identify modes other than those in Fig fTOh as 



"nonlinear modes" . 

We compare results at different variances with our discussion in £12.41 Some features of variance depen- 
dence are shared with the toy model in W2A\ the eigenvalue spectrum drifts and the corresponding modes 
rotate among themselves. However, the Abbott-Kepler modes also exhibit more complicated behavior. 
Fig [TT] shows projections of spike triggered stimuli onto two distinguished eigenvectors vx,+ and correspond- 
ing threshold crossing points. At low variance, the system crosses mostly one side of the threshold while 
the projections trace out the curvature of the threshold segment. As the variance increases, some crossing 
points begin to appear on the left side of the threshold. However, this does not overcome the expansion 
of a crossing point distribution on the right side, and the modes corresponding to this direction dominate. 
At high variance, there are many crossing points on both sides, reflected in the bimodal distribution of the 
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projections, previously seen in the toy models. Thus this is another example of how the results of reverse 
correlation analysis are affected both by the filter properties and by the interaction of the stimulus ensemble 
with the threshold geometry. 

So far, we have discussed only two dimensional dynamical models. We began with some artificial toy 
models with purely linear subthreshold dynamics, and proceeded to the minimal FitzHugh-Nagumo spiking 
model. We applied the lessons learned there to the more nonlinear Abbott-Kepler model. We will conclude 
with a partial analysis of the higher dimensional Hodgkin-Huxley model. 



6 Hodgkin-Huxley model 

Higher dimensional systems require nontrivial extensions of the methodology we have used with two dimen- 
sional systems. First, it is much harder to use the phase portrait to find a dynamical threshold, which could 
now be a multi-dimensional hypersurface rather than a curve. If we do not align the spike triggered stimuli 
according to the dynamical threshold, the obtained filters may include a distribution of time delays. For 
example, the reverse correlation analysis may result in broadened filters, 



/ M (t)= / drf^p^t-r), (35) 



where Pair) is a point spread function, depending on the input variance a 2 , as in Fig[ 

However, we can still gain some insights from / M . If p(t) has a narrow support, say less than a millisecond, 
it has negligible effects for larger time windows. It is also possible that some consequences of our analysis 
will still apply to / M with only a few reasonable assumptions. For example, if p a is of limited temporal 
extent, then the derivatives of / M can be simply 

fl n) ™ [dTf^(r) P(7 (t-T). (36) 



Therefore, as we discussed in <j2l the linear modes among {f^} should still be (approximately) closed under 
time differentiation. 

To demonstrate this, we use the four dimensional Hodgkin-Huxley model, and the following strategy: 
we select the linear modes out of the significant filters from covariance analysis. Now, since the STA is 
approximately their linear combination, the time derivatives of the STA should be written in terms of the 
linear combinations as well if they are closed under time differentiation. Fig [12] shows that this holds. In 
the low variance case, FigfTSa. the three linear modes chosen provide good fits to the time derivatives of the 
STA. The high variance case, Fig[T2"b. is affected by multi-spike effects and a filtering artifact, but it is also 
well fitted. Therefore, we can conclude that the time derivatives of the STA span the same space as three 
linear modes. One might expect a fourth linear mode due to the dimensionality of the model, but this is 
not as significant as the others; this agrees with previous covariance analysis that the HH model can be well 
described as a quasi-three dimensional system Aguera y Areas et al, 2003] . 



We note that it is not clear whether the linear modes in the two cases span the same feature space. This 
is difficult to ascertain since the point spread function can in principle depend on the stimulus variance. 
In Fig [T2"b . we see that while v\ of the high variance case can be fitted by the linear modes of the low 
variance, the other modes show a small deviation even around 5ms. This might indicate an interesting 
variance dependence as in ^2A\ but we will not pursue this issue in this paper. 



7 Summary 

Here we have investigated the meaning or interpretation of the features derived from the spike-triggered 
covariance method applied to two simple neuron models, the FitzHugh-Nagumo model, the minimal spiking 

7 In fact, it can be much more complicated than this. Since the time delay depends on the location of threshold crossing in a 
phase space, additional temporal correlations are introduced. Therefore, each filter may have a different point spread function, 
or filters may depend upon the time delay structure in a complicated way. We will discuss here only the simplest case. 
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neuron model, and the Abbott-Kepler model, a more faithful two-dimensional reduction of the Hodgkin- 
Huxley system. The power of white noise analysis is that it provides a data-driven method to reduce a 
high dimensional dynamical system to a functional model which captures the essential computation of the 
system in terms familiar to systems neuroscience: a receptive field which filters the stimulus, and a threshold 
function over the filtered stimulus. Our goal here was to analyze the output of a white noise analysis in 
terms of what it can reveal and how it depends upon the underlying dynamical system. In this, our approach 
is distinct from the elegant work of |Huys et al., 2006] where responses to white noise stimuli are used to fit 
the parameters of a conductance-based model. 

Dealing with simple two-dimensional systems, our observations of spiking dynamics in the phase plane 
motivated the following reduced model: dynamics in the phase plane are approximated by the perturbative 
expansion, in particular the linear approximation, and the system's nonlinearity is captured by a spiking 
threshold, determined for zero input, that extends through multiple dimensions. This model is a generaliza- 
tion of a basic filter-and-fire model to multiple dimensions, and extends it in two significant ways. One is 
the identification of the filters with the dynamical variables of the system; the other is the generalization of 
the concept of the threshold. 

The simplified model with linear dynamics and a multidimensional curved threshold was treated both 
analytically and phenomenologically, using numerical simulation and reverse correlation (covariance) analysis 
of spike triggered stimuli. This led to several insights, some of which are related to existing observations. 
First, for this case the feature space derived from covariance analysis is spanned by the linear kernels. Since 
the set of kernels is closed under time differentiation, the same feature space can be also spanned by a generic 
linear combination of the kernels, such as the spike triggered average and its n time derivatives. However, not 
every kernel contributes as a relevant feature. The threshold structure plays a role of selecting the relevant 
ones from the kernels: for example, a linear threshold selects a single filter and its time derivative. In general, 
a threshold spanning a d dimensional subspace will select d + 1 features from among linear combinations of 
the kernels. We show further for this model that the distribution of spike triggered stimuli in the covariance 
feature space corresponds to that of threshold crossing points, up to a suitable linear transformation. 

We also used this model to illustrate the effects of the interaction of the threshold nonlinearity with the 
stimulus ensemble. We show that a complex threshold geometry leads to a nontrivial variance dependence 
of the eigenmodes and eigenvalues of covariance analysis, and even more so, the spike-triggered average. 
In the linearized subthreshold model, the subspace of eigenmodes is not changed, but the spike-triggering 
ensemble may rotate through this subspace, leading to a variance-dependent spike-triggered average. This 
is one example of stimulus variance dependence in a nonlinear, non-adapting system |Yu and Lee, 2003} 



Borst et al., 2005 . We show this effect in the analysis of the Abbott-Kepler model neuron. 

The identification of a curved or dynamical threshold raises an issue in reverse correlation analysis regard- 
ing the determination of spike time. For the FitzHugh-Nagumo and Abbott-Kepler models, we compared 
two cases: when spike times are identified using a threshold in voltage, and when the spike times are found 
using the crossing of an identified curve that leads to spiking for zero current input. Surprisingly, the spike 
triggered averages showed remarkable differences. For a voltage threshold, the time delay from the earlier 
threshold crossing point blurs the estimated filters. It is shown that using the dynamical threshold leads to 
a better fitting of the estimated filters by the system's linear kernels, as well as improving the sharpness and 
substantially altering the shape of the STA. 

In all of our discussion here, we have concentrated on the subthreshold dynamics bringing the system to 
the point of threshold. We do not analyze, and our simplified models do not accommodate, the dynamics 
of the system immediately after a spike. In the phase plane picture, spiking reinjects the system into the 
subthreshold regime in a nonrandom way, affecting the subsequent probability flux to threshold, analogous to 
the one-dimensional reset of the integrate-and-fire neuron treated in |Paninski et al., 2003| . In our approach, 
we have assumed that the system has remained below threshold for long enough that its location in the sub- 
threshold space has been randomized by the driving current. This corresponds to an analysis of isolated spikes 
only, a simplification we have used before |Aguera y Areas et al., 20*03] and employed here for the covariance 
analysis. A number of works have treated this issue explicitly with a variety of methods: treating the inter- 
spike interval as the primary symbol |de Ruyter van S tcveninck and Bialck, 1988, Ricke 'eFal., 1997| , solv- 
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ing for the interspike interaction using the interspike intervals [Pillow and Si monccll i, 2003] , simultaneously 
solving for the linear kernel over stimulus history and spike history using the autoregressive moving average 
|Powers et al., 2005} Truccolo et al., 2005 and fitting parameters of an explicit model for the effective post- 
spike current [Kistler et al., 1997] |Keat et al., 2001[ |Paninski et al., 2004] |Pillow et al., 2005| . While the 
simplification we have used here allows us to find direct connections between the covariance modes without 
the confound of the interspike interaction [Agiiera y Areas et al. , 20 03 , Agii era y Areas and Fairhall, 2003] , 
it is clearly not a complete model for spiking responses. A first step toward a more complete spiking model 
may be to consider the perturbed subthreshold distributions induced by the influx of trajectories following 
spikes. This is in effect a mean field approximation, taking into account the overall spike rate for a given 
stimulus ensemble. Further steps could be taken by introducing a return map detcrministically relating 
the point of threshold crossing to a point of return into the subthreshold domain. The integrate-and-firc 
model is the most trivial implementation of such a return map; an equivalent map in multiple dimensions 
would reintroduce all spike trajectories into the subthreshold domain at an identified point (this is V = 
for integrate-and-fire) . Such a many-to-one map implies that the neuron's state is completely reset by a 
spike, which is incorrect for neurons with slow conductances that modulate spike afterpotentials. An less 
degenerate map seems more appropriate for such cases. Another important step is the addition of these slow 
conductances. Using our formalism, such conductances may be representable simply as additional dimensions 
of threshold curvature with corresponding longer-timescale stimulus filters. 

White noise analysis allows the derivation of intuitive functional models for neural computation: what 
does the neuron compute? In this paper we have drawn concrete correspondences between the components 
of these functional models and parameters of the underlying dynamical system. 



A Volterra expansion of a dynamical system 

Let's consider an n-dimensional dynamical system, perturbed by an external input I(t) as Eq (|15p . For 
convenience, we introduce a dimensionless expansion parameter rj, with which we will expand the equation. 
More precisely, Eq (fT5|) becomes 



dtVk = fk{yi,V2, ■••,yjv) +r)5 k0 I(t), (37) 
We also have the perturbation expansion of yj~ as 

» fc = yJP ) + rf ) + 7 a yf ) + -. (38) 

By plugging Eq ((38]) into Eq 137]), we obtain 

« ( W _l 2 (2) , \ Tf .s x . %(?/ (0) ) (1) 

dt [vy k + v v k + ■ ■ ■ ) = m{t)ho + — ww 

i d 2 A-0/ 0) ) 2 (1) (1) , df k (y(°y) 2 (2) 

By comparing two sides order by order, we obtain a series of equations satisfied by the perturbative expansion 
at each order as 



T^kmV^m — I^kO, T^km = Skmdt 



which is Eq (fTB]) . 

t>. y(2) - 1 Hl w (i) u (i) Hi - d 2 -^ (0) ) (391 

^kmilm 2! kmnilm tin 1 MJ -kmn > K ^) 

which is Eq ([23]) . and so on. 

These equations can be solved recursively by using a kernel e^ 1 ^ = V^ 1 which therefore satisfies 
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Now the solution of the first order can be written as 

y k 1] (t) = J dt'e^it-t'W). 

where = e^. Also from Eq ([39]) . 



(2) 

Vk 



I J dt'ei)\t-t') (H lrnn J d Sl ds 2 eV{t' - s x )$\t' - s 2 )I{si)I{s 2 ) 

= 2~! j ds T- ds 2 ( f ~ Slit- S 2 )I(si)I(s 2 ), 

e[ 2 \suS2) = J dt'H lmn e$(t')e$( Sl - t')e^(s 2 - t% 

which is Eq (|24p . Since this procedure can be carried out to higher orders, the higher order kernels are outer 
products of the linear kernels. 

In addition, the first order kernels are related to each other by simple time differentiation. To show this, 
let's rewrite Eq (TTg]) as follows 

e$(t-t')=V-l5 n0 5(t-t'). (40) 
Now the derivatives can be computed as 

d t e$ = d t -D-'j n0 6(t-t') 

= (d t S m i - J ml + J m i)T>Y^5 n0 5(t - t') 

= V m iV-^5 n0 5{t - t') + J ml e\ 1] 

= 6 m0 S{t - t') + J m/ e ; (1) (41) 

This shows that the non-singular part of the first order kernels can be obtained by the linear combination 
of the derivatives of other kernels. 

B A linear model from the FitzHugh-Nagumo system 

Here we derive the linearization of the FitzHugh-Nagumo model, beginning with Eqs. ([5]) and (J6j> . The fixed 
point can be obtained by simultaneously solving the nullcline equations. We denote it by (Vb, Wo), and we 
expand the system around this point. 
First, the Jacobian is 

F'{V )/^ -l/V 
1 -b 



J = 



where F(V) — V(l — V)(a + V), and this defines a linear system 

F'(Vn) 1 

d t V = -^-V--W + I(t) 

1p 1p 



d t W = V -bW. 



J has eigenvalues A±, 



A± = - 



f-(V )± yJf+(V Q ) 2 -lP 



-b+ [/+(7o)±V/+W-V 

-b + k±. 



(42) 
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where f±(V) = (F'(V) ± bijj)/2. As in Eq. (gD]), J is diagonalized by a matrix U, 

U= 1 ( 1 X ~ +b ) 
\+-X-\-l -(A+ + 6) ; ' 

From Eq. (|22p . we obtain the first order kernels which solve the linear system 

e$\t) = e- bt d t S(t)H(t), 



(43) 



e«(t) = e- 6t 5(f)ff(t) ) 



where S(t) = (e K +* — e K ~ t )/{n + — «;_). with our choice of parameters is drawn in Fig[2k- 

C Derivation of the Abbott-Kepler model 

In this section, we show the derivation of a two-dimensional neuron model used in Sj5] For further details, 
we refer to the original paper |Abbott an d Kcplc r~1990| . 

We begin with a Hodgkin-Huxley equation, which is defined by Eq. (fTJ) and the following parameters: 

9l = 9l, 9k = 5ifn 4 , g Na = g Na m 3 h. (44) 

T z(V)^ = z(V) - z, t z = \— , z= — z = m,n,h. (45) 

at a z + p z a z + (3 Z 

1(V + AO) 

p m = 4exp[.0556(V + 65)], 



l-exp[-.l(V + 40)]' 
cp[.05(V + f 
.01(V + 55) 



a h = .07«p[.05(V + 65)] ) & = t + 35)] , (46) 

@ n = .125exp[.0125(V + 65)]. 



l-exp[-.l(V + 55)]' 

Now, the key observation is that r m is much smaller than r„^ while r/j and r„ are mutually comparable. 
Therefore, m can be approximated by its value at equilibrium, m(V), and h and n can be represented by 
the same equilibrium voltage, U. In other words, 

mram(V), h^h(U), n^n(U). 

Using this, we obtain Eq. ([33]) as 

i=L,K,Na 

« 9l(V - E L ) + g K n(U) 4 (V - E K ) + g Na m(V) 3 h(U)(V - E Na ) 
= -.f(V,U). 

An equation for U, Eq. (|34[) . is obtained by requiring that time dependence of the active current F due to 
h and n in the Hodgkin-Huxley model is mimicked by h(U) and n(U). This implies 

dF dh dF dn fdf dh df dn \ dU 



, (47) 
dh dt dn dt V dh dU dn dU J dt K ' 



Again, we approximate Eq. (|45[) in the same way, 

dz 1 



(z(V)-z(U)), z = h,: 



dt (V) 

Plugging this in Eq. 1(17)). we can solve for dU/dt as a function of V and [/, which we denoted by g(V, U) in 
Eq. (J34J) as 

ffjVa(U - E Na )m{Vf(h{V) - h(U))/r h (V) + 4g K (V - E K )n{Uf{n{V) - h(U))/r n (V) 



g(v,u) 



g Na (V - E Na )m{yfh'{U) + Ag K {V - E K )n(Ufn'(U) 
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Figure 1: (a) Phase plane of the FitzHugh-Nagumo neuron model with trajectories and threshold (dashed). 
Parameters are a — 0.1, b — 0.5, c = —0.4 and ip — 0.015. (b) A single extended trajectory on the FN 
phase plane driven by white noise input, showing multiple spikes. The threshold (dashed) still makes an 
approximate boundary of the spiking and non-spiking portions of the trajectory (gray). 
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Figure 2: (a) First order kernels of the FN neuron model, drawn as niters, i. e. in — t.(b) Various thresholds 
for the linear system. V (vertical), W (horizontal), linear (dashed), piecewise linear (dotted), and curved 
(solid curve) thresholds are used, (c) Spectra of covariance matrices for the thresholds. The shaded box 
represents the level of error from finite sample size and dimension, (d) STAs and covariance modes of the 
linear system with various thresholds. The non-singular part of each mode is fitted using least squares to 
a linear combination of the (orthogonalized) first order kernels as Vfu — cyS^ + cw^y) ■ The kernels are 
normalized excluding the ^-function component. The coefficients (cy, cw) are displayed above each plot and 
the gray line is a fitted function. For each case we show the projection of spike-triggered current histories 
onto either the leading two negative modes, v± and V2, or v\ and i>+, depending on the threshold. 
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Figure 3: (a) A diagram of a single filter model. The disc of radius a represents the prior Gaussian 
distribution. The dotted line is the model threshold imposed at zq = 8f. The thick line is a distribution of 
spike triggered stimuli, and the white dot represents its average, the STA. (b) A diagram of a model whose 
threshold is the union of two lines, each imposed at a distance 9f, g from the center. Two gray dots denote 
the STAs for each segment. 
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Figure 4: Two examples of variance dependence in the linear system model, (a) Spectra of the covariance 
matrix for several different variance stimuli, (b) Significant eigenmodes. (c) Projection of spike-triggering 
stimuli onto the filter subspace defined by for the linear and v\^+ for the piecewise linear threshold 
respectively, colored according to the stimulus variance as in the legend of (b). The circled points are 
the centers of the distributions, corresponding to the respective STAs, colored according to the respective 
distribution. 
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Figure 5: Effect of a dynamical threshold in the FN model, (a) STA of the full FN system with a constant 
V threshold. As before, a linear fit by the kernels is also shown, (b) Comparison of spike triggered stimuli 
with different threshold choices in the FN model, (c) Temporal point spread function induced by the choice 
of threshold. At is the time delay from crossing the dynamical threshold to the fixed V threshold, (d) Time 
delay of each spike triggered stimulus plotted against the value of W at its threshold crossing point. 
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Figure 6: (a) STAs of a FN model neuron, and the first (gray) and second (dashed) order approximation. (b) 
Spectra of covariance matrices of the FN model and approximations, (c) STAs and covariance modes for the 
first and second order approximations and the full FN system. 
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Figure 7: (a) Spike triggered stimuli projected into a feature space defined by v\.+ for the first order case 
and full system, and for the second order system, (b) The density of threshold crossing points in the 
first order (solid), second order (dashed), and full (dotted) systems, plotted along the threshold curves in 
the V — W plane. 
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Figure 8: (a) Phase plane of the Abbott-Kepler model with trajectories, nullclines, and a threshold, in the 
absence of input, (b) AK phase plane with the injected noisy input. 
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Figure 9: (a) Spike triggering current histories with two different threshold schemes in the AK model, (b) 
STAs of the AK model with the fixed V = — 40mV threshold (dashed) and dynamical threshold (solid). 
Also, they are least-squares fitted by linear combinations of ev,u (gray)- 
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Figure 10: (a) Normalized first order kernels ev,u of the AK model and the component of ejj orthogonal 
to ey, normalized (gray), (b) Eigenvalue spectrum of the covariance matrix of the AK model with each 
variance, oq = 13pA. (c) Eigenmodes of the covariance matrices. As previously, «i 2 and i> + correspond to 
the two smallest negative eigenvalues and the largest positive one. 
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Figure 11: Projection of spike triggering stimuli onto V\ and v + in the AK model. Below are threshold cross- 
ing points in the phase plane of each case, marked by crosses, V < — 67mV, and circles, when V > — 67mV, 
over the dynamical threshold (dashed). The other gray curves are the nullclines, with their intersection, the 
fixed point, marked by a circle. 
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Figure 12: Covariance bases of the linear feature space of the Hodgkin-Huxley model recovered at (a) low 
and (b) high variance. The three most significant eigenvalues are denoted with the same colors as the 
corresponding linear eigenmodes. (c) The STA and linear modes of the high variance case and their fits 
using the linear modes of the low variance case in the HH model. 
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